*********************************************************
*********************************************************
*                  RELIGIOSITY PROJECT					*
*					 PROGRAMS DOFILE                    *		
***      created by Hao Wang on Aug 27 Aug 2023      ****
** last & once for all modified by Hao on Aug 31 2023  **
*********************************************************
*********************************************************


***********
**# Index
***********

* This do file defines helpful programs to produce the correlation tables and regression tables

	* Programs:

		* correlations:
			* Produces correlation Tables in the paper
			
		* regressions:
			* Produces OLS regression tables in the paper
			
		* mlogitreg:
			* Produces Multinomial Logit regression tables in the paper
			

* Correlation Table
cap program drop	correlations
program define 		correlations
		
	cap erase "$tables/$tablename.tex"	
	
	foreach outcome in  $outcomes_corr {
		eststo: estpost corr `outcome' $covariates_corr
	}	
	

	esttab using "$tables/$tablename.tex", replace ///
			b(%5.3f) se(%5.3f) nogaps star(* 0.1 ** 0.05 *** 0.01) ///
			nomtitles nocons nolz title("Correlation Matrix 3 (KLPS 3)" ) ///
			label se noobs nonumbers collabels(none) longtable fragment nomtitles booktabs ///
			mgroups("\vspace{-.25cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{Panel A: Sample - $panelA}}} \vspace{.2cm}\\ \mainHeader$header", span) 
										
	estimates clear	
	
	
	loc heterogeneity $het1 $het2 
	foreach het of local heterogeneity {

		if     "`het'" == "$het1" {
				local letter "B" 
				local label "Panel B: Sample - $panelB"
				local table age
		}
		else if "`het'" == "$het2" {
				local letter "C"
				local label "Panel C: Sample - $panelC"
				local table age
		}
				
		foreach outcome in   $outcomes_corr {
		eststo: estpost corr `outcome' $covariates_corr if `het'==1 
		
		}		
	
	esttab using "$tables/$tablename.tex", append ///
			b(%5.3f) se(%5.3f) nogaps star(* 0.1 ** 0.05 *** 0.01) ///
			nomtitles nocons nolz title("Correlation Matrix 3 (KLPS 3)" ) ///
			label se noobs nonumbers collabels(none) longtable fragment nomtitles booktabs ///
			mgroups("\vspace{-.25cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{`label'}}} \vspace{.2cm}\\ \mainHeader$header", span) 
											
	estimates clear
				
	}	
		
end
		

		
* OLS Regression Table
cap program drop 	regressions
program define 		regressions
		
		
	cap erase "$tables/$tablename.tex"
	estimates clear
	
	foreach outcome of varlist $outcomes {		
		
		* save percentiles for trimming					
		sum `outcome', d
		
		*control mean & sd - not trimmed
		sum `outcome' if psdp_treat == 0 [aw=con_weight] 
		loc controlMean = r(mean)
		loc controlSD = r(sd)
				
		* Pooled OLS estimator
		eststo: reg `outcome' psdp_treat $controls  [pw=con_weight] ,  cluster(con_psdpsch98) 
		loc treat = _b[psdp_treat]
		estadd scalar control_mean = `controlMean'
		estadd scalar control_SD   = `controlSD'
		estadd loc Controls Yes, replace
		loc teffect = 100*(`treat'/`controlMean')
		estadd scalar t_effect = round(`teffect', .01)	
				
		* saving number of unique pupid observations in regression sample
		preserve
		keep if e(sample)
		bysort pupid (pupid): gen c`outcome' = 1 if _n==1 
		summ c`outcome'
		estadd scalar N_ind = r(N)
		restore	
	}
					

				

		esttab using "$tables/$tablename.tex", ///
				b(%5.3f) se(%5.3f) r2 nogaps star(* 0.1 ** 0.05 *** 0.01) ///
				nomtitles nocons nolz label se noobs nonumbers collabels(none) replace longtable fragment nomtitles booktabs ///
				mgroups("\vspace{-.4cm}\\ \toprule \vspace{-.4cm} \\&\multicolumn{7}{c}{\emph{\textbf{Panel A: Sample - $panelA}}} \\ \toprule   \\ \mainHeader$header", span) ///
				keep(psdp_treat) stats(control_mean control_SD t_effect r2 N_ind N, ///
				labels("Control Mean" "Control SD" "Treatment Effect (\%)" "R-squared" "Number Individuals" "Number Observations") fmt(2 2 2 2 0)) ///
				substitute("Standard errors in parentheses" " ""\sym{*} \(p<.10\), \sym{**} \(p<.05\), \sym{***} \(p<.01\)" " ") ///					

		estimates clear
			
		
		loc heterogeneity 	$het1 $het2 
		foreach het of local heterogeneity {

						if     "`het'" == "$het1" {
								local letter "B" 
								local label "Panel B: Sample - $panelB"
								local table age
						}
						else if "`het'" == "$het2" {
								local letter "C"
								local label "Panel C: Sample - $panelC"
								local table age
						}
				
			foreach outcome of varlist $outcomes {		
						su `outcome', d
						*control mean & sd - not trimmed
						sum `outcome' if psdp_treat == 0 & `het' == 1 [aw=con_weight] 
						loc controlMean = r(mean)
						loc controlSD = r(sd)
					* Pooled OLS estimator by gender
						eststo: reg `outcome' psdp_treat $controls  if `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98)  
							loc treat = _b[psdp_treat]
						estadd scalar control_mean = `controlMean'
						estadd scalar control_SD = `controlSD'
						estadd loc Controls Yes, replace
						loc teffect = 100*(`treat'/`controlMean')
						estadd scalar t_effect = round(`teffect', .01)	
						** saving number of unique pupid observations in regression sample
							preserve
							keep if e(sample)
							bysort pupid (pupid): gen c`outcome' = 1 if _n==1 
							summ c`outcome'
							estadd scalar N_ind = r(N)
							restore		
					}
						
		# delimit ;
				esttab using "$tables/$tablename.tex", 
				b(%5.3f) se(%5.3f) r2 nogaps star(* 0.1 ** 0.05 *** 0.01)
				nomtitles nocons nolz 
				label se noobs nonumbers collabels(none) append longtable fragment nomtitles booktabs
				mgroups("\vspace{-.25cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{`label'}}} \vspace{-.2cm}\\", span) 
				keep(psdp_treat) 
				stats(control_mean control_SD t_effect r2 N_ind N, 
						labels("Control Mean" "Control SD" "Treatment Effect (\%)" "R-squared" "Number Individuals" "Number Observations") fmt(2 2 2 2 0))
				substitute("Standard errors in parentheses" " "
										"\sym{*} \(p<.10\), \sym{**} \(p<.05\), \sym{***} \(p<.01\)" " ")						
			;
			#delimit cr		
			estimates clear
			eststo clear	

		}	 
		
		end		
		
		
* Multinomial Logit Regression Table
cap program drop 	mlogitreg
program define 		mlogitreg
		
		
	cap erase "$tables/$tablename.tex"
	estimates clear
	
		foreach outcome of varlist $outcomes {		
				* save percentiles for trimming					
					sum `outcome', d
				*control mean & sd - not trimmed
					sum `outcome' if psdp_treat == 0 [aw=con_weight] 
					loc controlMean = r(mean)
					loc controlSD = r(sd)
				* Pooled OLS estimator
					eststo: mlogit `outcome' psdp_treat $controls  [pw=con_weight] ,  cluster(con_psdpsch98) 
					  loc treat1 = [#2]_b[psdp_treat]
					  loc treat2 = [#3]_b[psdp_treat]
					  estadd scalar rsq = e(r2_p)
						estadd scalar   control_mean = `controlMean'
						estadd scalar   control_SD   = `controlSD'
						estadd loc 		Controls Yes, replace
							loc teffect = 100*(`treat1'/`controlMean')
						estadd scalar t_effect1 = round(`teffect', .01)	
							loc teffect = 100*(`treat2'/`controlMean')
						estadd scalar t_effect2 = round(`teffect', .01)	
						** saving number of unique pupid observations in regression sample
						preserve
						keep if e(sample)
						bysort pupid (pupid): gen c`outcome' = 1 if _n==1 
						summ c`outcome'
						estadd scalar N_ind = r(N)
						restore	
						}
				
		# delimit ;
			esttab _all using "$tables/$tablename.tex", 
			noomitted modelwidth(13) unstack compress
			b(%5.3f) se(%5.3f) r2 nogaps star(* 0.1 ** 0.05 *** 0.01)
			nomtitles nocons nolz 
			label se noobs nonumbers collabels(none) replace longtable fragment nomtitles booktabs
			mgroups("\vspace{-.4cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{Panel A: Sample - $panelA}}} \\ \toprule   \\ \mainHeader$header", span) 
			keep(psdp_treat) 
			stats(control_mean control_SD t_effect1 t_effect2 rsq N_ind N, 
					labels("Control Mean" "Control SD" "Treatment Effect 1 (\%)"  "Treatment Effect 2 (\%)" "R-squared" "Number Individuals" "Number Observations") fmt(2 2 2 2 2 0))
			substitute("Standard errors in parentheses" " "
									"\sym{*} \(p<.10\), \sym{**} \(p<.05\), \sym{***} \(p<.01\)" " ")						
		;
		#delimit cr	
		estimates clear
			
		
		loc heterogeneity 	$het1 $het2 
		foreach het of local heterogeneity {

						if     "`het'" == "$het1" {
								local letter "B" 
								local label "Panel B: Sample - $panelB"
								local table age
						}
						else if "`het'" == "$het2" {
								local letter "C"
								local label "Panel C: Sample - $panelC"
								local table age
						}
				
			foreach outcome of varlist $outcomes {		
						su `outcome', d
						*control mean & sd - not trimmed
						sum `outcome' if psdp_treat == 0 & `het' == 1 [aw=con_weight] 
						loc controlMean = r(mean)
						loc controlSD = r(sd)
					* Pooled OLS estimator by gender
						eststo: mlogit `outcome' psdp_treat $controls  if `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98)  
						 loc treat1 = [#2]_b[psdp_treat]
						 loc treat2 = [#3]_b[psdp_treat]
							estadd scalar   control_mean = `controlMean'
							estadd scalar   control_SD   = `controlSD'
							estadd loc 		Controls Yes, replace
								loc teffect = 100*(`treat1'/`controlMean')
							estadd scalar t_effect1 = round(`teffect', .01)	
								loc teffect = 100*(`treat2'/`controlMean')
							estadd scalar t_effect2 = round(`teffect', .01)	
						** saving number of unique pupid observations in regression sample
							preserve
							keep if e(sample)
							bysort pupid (pupid): gen c`outcome' = 1 if _n==1 
							summ c`outcome'
							estadd scalar N_ind = r(N)
							restore		
					}
						
		# delimit ;
				esttab using "$tables/$tablename.tex", 
				noomitted modelwidth(13) unstack compress
				b(%5.3f) se(%5.3f) r2 nogaps star(* 0.1 ** 0.05 *** 0.01)
				nomtitles nocons nolz 
				label se noobs nonumbers collabels(none) append longtable fragment nomtitles booktabs
				mgroups("\vspace{-.25cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{`label'}}} \vspace{-.2cm}\\", span) 
				keep(psdp_treat) 
				stats(control_mean control_SD t_effect1 t_effect2 r2_p N_ind N, 
						labels("Control Mean" "Control SD" "Treatment Effect 1 (\%)" "Treatment Effect 2 (\%)" "R-squared" "Number Individuals" "Number Observations") fmt(2 2  2 2 2 0))
				substitute("Standard errors in parentheses" " "
										"\sym{*} \(p<.10\), \sym{**} \(p<.05\), \sym{***} \(p<.01\)" " ")						
			;
			#delimit cr		
			estimates clear
			eststo clear	

		}	 
end		
